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Abstract. The aim of this paper is to extend the global error estimation and control 
addressed in Lang and Verwer [SIAM J. Sci. Comput. 29, 2007] for initial value prob- 
lems to finite difference solutions of parabolic partial differential equations. The classical 
ODE approach based on the first variational equation is combined with an estimation of 
the PDE spatial truncation error to estimate the overall error in the computed solution. 
Control in a discrete L2-norm is achieved through tolerance proportionality and mesh re- 
finement. Numerical examples are used to illustrate the reliability of the estimation and 
control strategies. 



1. Introduction 

We consider initial boundary value problems of parabolic type, which can be written as 

(1) d t u(t,x) = f(t,x,u(t,x),d x u(t,x),d xx u(t,x)) , t G (0, T] , x G Q C R d , 
equipped with an appropriate system of boundary conditions and with the initial condition 

(2) u(0, x) = u (x) , x G ft. 

The PDE is assumed to be well posed and to have a unique continuous solution u(t, x) which 
has sufficient regularity. 

The method of lines is used to solve (pQ) numerically. We first discretize the PDE in space 
by means of finite differences on a (possibly non-uniform) spatial mesh Qh and solve the 
resulting system of ODEs using existing time integrators. For simplicity, we shall assume 
that this system of time- dependent ODEs can be written in the general form 

U' h (t) = F h (t,U h (t)), tG(0,T], 

u h (o) = [4,0, 

with a unique solution vector Uh(t) being a grid function on Q h . Let 

(4) R h : u(t, ■ ) -> R h u(t) 

be the usual restriction operator defined by Rh.u(t) = (u(t, x{), ■ ■ ■ , u(t, xn)) t , where Xi G flh 
and N is the number of all mesh points. Then we take as initial condition Uh,o = Rhu(0). 

To solve the initial value problem ([3]), we apply a numerical integration method at a certain 
time grid 

(5) = t < tx < ■ ■ ■ < t n < ■ ■ ■ < £ M _i < t M = T , 
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using local control of accuracy. This yields approximations Vh(t n ) to Uh(t n ), which may be 
calculated for other values of t by using a suitable interpolation method provided by the 
integrator. The global time error is then defined by 

(6) e h (t) = V h (t) - U h (t) . 

Numerical experiments in [lj for ODE systems have shown that classical global error estima- 
tion based on the first variational equation is remarkably reliable. In addition, having the 
property of tolerance proportionality, that is, there exists a linear relationship between the 
global time error and the local accuracy tolerance, eh(t) can be successfully controlled by a 
second run with an adjusted local tolerance. Numerous techniques to estimate global errors 
are described in [8]. 

In order for the method of lines to be used efficiently, it is necessary to take also into 
account the spatial discretization error. Defining the spatial discretization error by 

(7) Vh(t) = U h (t) - R h u(t) , 

the vector of overall global errors Eh(t) = Vh(t) — Rhu(t) may be written as sum of the global 
time and spatial error, that is, 

(8) E h {t) = e h {t) + Vh {t) . 

It is the purpose of this paper to present a new error control strategy for the global errors 
Eh{t)- We will mainly focus on reliability. So our aim is to provide error estimates Eh(t) « 
Eh(t) which are not only asymptotically exact, but also work reliably for moderate tolerances, 
that is for relatively coarse discretizations. 

The global errors are measured in discrete L2-norms. A priori bounds for the global error 
in such norms are well known, see e.g. [5], [9] . However, reliable a posteriori error estimation 
and efficient control of the accuracy of the solution numerically computed to an imposed 
tolerance level are still challenging. We achieve global error control by iteratively improving 
the temporal and spatial discretizations according to estimates of e^(t) and rjh(t). The global 
time error is estimated and controlled along the way fully described in [3]. To estimate the 
global spatial error, we follow an approach proposed in pQ (see also [B]) and use Richardson 
extrapolation to set up a linearised error transport equation. 

2. Spatial and time error 
By making use of the restriction operator Rh, the spatial truncation error is defined by 

(9) a h (t) = (R h u)'(t) - F h (t, R h u(t)) . 

From ([3]) and , it follows that the global spatial error rjh (t) representing the accumulation 
of the spatial discretization error is the solution of the initial value problem 

V' h (t) = F h (t,U h {t))-F h {t,R h u{t))-a h (t), te(0,T], 

ifo(O) = o. 

Assuming Fh to be continuously differentiate, the mean value theorem for vector functions 
yields 

V' h (t) = du.F^Uhity^-a^ + OMt) 2 ), te(0,T], 
Vh(0) = 0. 



Global Error Estimation and Control 3 

With Vhit) being the continuous extension of the numerical approximation to ([3]), the residual 
time error is defined by 

(12) r h (t) = V^t)-F h (t,V h (t)). 

Thus the global time error eh{t) fulfills the initial value problem 

e' h (t) = F h (t, V h (t)) - F H (t, U h (t)) + r h (t) , t 6 (0, T] , 
c fc (0) = 0. 
Again, the mean value theorem yields 

e' h (t) = d Uh F h (t,V h (t))e h (t)+r h (t) + 0(e h (t) 2 ), t e (0,T\, 
e h (0) = 0. 

Apparently, by implementing proper choices of the defects au{t) and rhit), solving (TiT]) and 
fTl4|) will in leading order provide approximations to the true global error. The issue of how 
to approximate the spatial truncation error and the residual time error will be discussed in 
the next sections. 

3. Estimation of the residual time error 

We assume that the time integration method used to approximate the general ODE system 
([3]) is of order p < 3. Following the approach proposed in jl] we define the interpolated 
solution Vhit) by piecewise cubic Hermite interpolation. Let Vh, n = Vh(t n ) and F^ n = 
Fh(t n , Vh >n ) for all n — 0, 1, . . . , M. Then at every subinterval [t n , £ n +i] we form 

(15) V h (t) = V h , n + A n (t-t n ) + B n (t-t n ) 2 + C n (t-t n ) 3 , t n <t<t n+1 , 
and choose the coefficients such that V^(t n ) = F h ^ n and V^(i n+ i) = F hn+ i. This gives 

(16) V h (t n + 9r n ) = v (9)V hin + Vi(9)V h;n+1 + r n w (9)F htn + T n w x {9)F Kn+1 
with < 9 < 1, r n = t n+ \ — t n , and 

(17) v (9) = (1 - 9) 2 {l + 29), Vl (9) = 9 2 (3 - 29), w (9) = (1 - 9) 2 9, w x {9) = 9 2 (9 - 1). 

Now let Yh(t) be the (sufficiently smooth) solution of the ODE (J3]) with initial value Y(t n ) = 
Vh, n - Then the local error of the time integration method at time t n+ i is given by 

(18) le n+1 = V h>n+1 - Y h (t n+1 ) = 0(t? +1 ). 
Combining ffTB"]) and ffTS]) gives 

V h (t n + 9r n ) - Y h (t n + 9r n ) = 

vi {9)le n+1 - Y h {t n + 9r n ) + v {9)Y h {t n ) + v 1 (9)Y h (t n+1 ) 

( 19) 

+T n w (9)F h (t n , Y h (t n )) + T nWl (9)F h (t n+1 , Y h (t n+1 )) 
+T n Wi(9)(F h (t n+1 , 14,n+i) - F h (t n+ i, Y h (t n+1 ))), 
and by Taylor expansion we obtain 

(20) V h (t n + 9r n ) - Y h (t n + 9r n ) = v 1 {9)le n+1 + ^(29 3 - 9 2 - 9 A )r A n Y h (A) (t n ) + 0(r p n +2 ) . 

Recalling Y^it) = F h {t, Y h (t)) for t G (t n , t n+ i] and rewriting the residual time error as 

(21) r h (t) = V{£t n + 9r n ) - Y^{t n + 9r n ) + F h (t, Y h (t)) - F h (t, V h (t)) , 
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with 9 = {t — t n )/r n , we find by differentiating the right hand side of (120]) 

(22) r h (t n + 6r n ) = 6(9 - 6 2 ) 1 -^ + ±-(36 2 - 6 - 2^)r n 3 F /l (4) (t„) + 0(r^) . 

T n 12 

Setting 6 = 1/2 in will reveal 

(23) r h (t n+1/2 ) = |^±i + ^(rr 1 ) . 

Thus the cubic Hermite defect halfway the step interval can be used to retrieve in leading 
order the local error of any one-step method of order 1 < p < 3 (see also [I], Section 2.2). 
Following the arguments given in |4J, Section 2.1, we consider instead of (j!4p the step size 
frozen version 

, OA s e' h (t) = du h F h (t n , V h , n ) e h (t) + %r h (t n+ i), t e (t n , t n+1 ), n = 0, . . . ,M-1, 

(24) e h (0) = 

to approximate the global time error Ch(t). Using 

1 T 

(25) Vh(t n+1 / 2 ) = -iVh,n + Vh,n+l) + -^{Fh,n ~ Fh, n +l) 

and 

3 1 

— {Vh,n+1 — Vh n ) — ~( 
It 4 

we can compute the residual time error halfway the step interval from (|12p 

fh{tn+l/2) = ^(Vh,n+1 — Vh, n ) ~ \(Fh, n + -^h,n+l) 



(26) V^(t n+ i/2) — — (Vh,n+1 ~ Vh,n) — ~(Fh,n + F^ n +l) 



(27) 



-F h U n +i> \ (Yh,n + Vh,n+l) + \ {Fh,n ~ Fh >n -i 



Remark 3.1. From 0221) we deduce 



1 r tn + 1 1p 
(28) - / r h (t) dt = ^ + 0(r^) 



showing, in the light of (j23|) . that ^rh(t n+1 / 2 ) is in leading order equal to the time-averaged 
residual. Thus, we can justify the use of the error equation f[24l) without the link to the first 
variational equation. 

4. Estimation of the spatial truncation error 

An efficient strategy to estimate the spatial truncation error by Richardson extrapolation 
is proposed in PQ. We will adopt this approach to our setting. 

Suppose we are given a second semi-discretization of the PDE system (CQ), now with dou- 
bled local mesh sizes 2h, 

U' 2h (t) = F 2h (t,U 2h (t)), 4 6(0,21, 

Mo) = u 2hfi . 

In practice, one first chooses Q 2 h and constructs then Qh through uniform refinement. The 
following two assumptions will be needed, (i) The solution U 2 h(t) to the discretized PDE 
on the coarse mesh Q 2 h exists and is unique, (ii) The spatial discretization error r)h(t) is of 
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order q with respect to h. We define the restriction operator R% h from the fine grid Qh to 
the coarse grid Q 2 h by the identity R 2 h = R 2 hRh and set 

(30) rjl(t) = R h 2h rj h (t), U c h (t) = R% h U h (t), V h c (t) = R% h V h (t) . 
From the second assumption it follows that 

(31) V c h (t) = 2-« V2h (t) + 0(h« +1 ) 
and therefore 

(32) R 2h u{t) = J-^U^(t) - i^U 2h {t) + O(h^) . 
The relation U%(t) — U2h{t) = Vh(t) ~ V2h(t) together with (T3T|) gives 

(33) um - u 2h {t) = l -^-m h {t) + o(h« +1 ) . 

The spatial truncation error on the coarse mesh Q 2 ^ is analogously defined to as 

(34) a 2h (t) = (R 2h u)'(t) - F 2h (t, R 2h u{t)) . 

Substituting R 2 hii(t) from (1321) into the right-hand side, using the ODE system (1291) to 
replace U 2h (t), and manipulating the expressions with (f3"3"|) we get 

(35) 

+^(F 2h (t,U 2h (t)- V2h (t) + 0(h<i +1 j) - F 2h (t,U 2h (t)j) + 
Taylor expansions yield 

(36) a 2h (t) = ^^((U c h y(t)-F 2h (t,U c h (t)))+0(h^). 

Analogously to ®, we set e%(t) = V£(t) - U c h (t). Substituting (E/£)'(t) by R% h F h (t,U h (t)) 
and using again Taylor expansion it follows that 

a 2h (t) = ^(R% h F h (t,V h (t))-F 2h (t,VZ(t)))+0(h^ 

(37) v 7 
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21- 



T (R h 2h (d Uh F h (t,V h (t))e h (t)) - d Uh F 2h (t^(t))el(t) ) + 0{e h {tf) . 



Assuming the term on the right-hand side involving the global time error to be sufficiently 
small, we can use 

(38) a 2h (t) = [R h 2h F h {t, V h (t)) - F 2h (t, V£(f))) 

as approximation for the spatial truncation error on the coarse mesh. To guarantee a suitable 
quality of the estimation (13"8l we shall first control the global time error for attempting that 
afterwards the overall error is dominated by the spatial truncation error (see Section [H]). 

An approximation 6th{t) of the spatial truncation error on the (original) fine mesh is 
obtained by interpolation respecting the order of accuracy (see Section [5]). Thus, to approx- 
imate the global spatial error i]h(t) we consider instead of (TTTj) the step-size frozen version 

fj' h (t) = du h F h (t n ,V ht n)fj h (t) -a h (t), te (t n ,t n+ i], n = 0,...,M-l, 

(o9) fj h (0) = 
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Remark 4.1. If an approximation e^(t) of the global time error has already been computed, 
we could make use of Uf t {t) ~ Vfi(t) + e c h (t) to obtain a better approximation of at2h(t) from 
( |36|) . However, we have found by experiments that even in the case when the global time 
error was not small, using the step size frozen equations (124"|) and (1391) together with (1361) to 
approximate the global time and spatial error does not yield a significantly better approx- 
imation. Since in practise the use of formula (}3"oT) requires additional function evaluations, 
equation (1551) appears to be more efficient. 

Remark 4.2. We note that special care has to be taken in the handling with the spatial 
truncation error at the boundary when derivative boundary conditions are present. This 
requests interpolation adopted to the correct order of accuracy, see [TJ. 

5. The example discretization formulas 

In order to keep the illustration as simple as possible we restrict ourselves to one space di- 
mension. For the spatial discretization of ([I]) we use standard second-order finite differences. 
Hence we have q = 2. The discrete L2-norm on a non-uniform mesh 

(40) x < xi < . . . < x N < xn + i , hi = Xi — , i — 1, . . . , N + 1 , 
for a vector y = . . . , i/n) T £ ^ N is defined through 

/ a-i \ II II 2 hi + hi+l 2 

( 41 ) \\v\\ =2^ — 2 — Vi ' 

i=i 

Here, the components y and Vn+i which are given by the boundary values are not considered. 

The example time integration formulas are taken from [1]. For the sake of completeness we 
shall give a short summary of the implementation used. To generate the time grid ([3]) we use 
as an example integrator the 3rd-order, A-stable Runge-Kutta-Rosenbrock scheme ROS3P, 
see [21 13] for more details. The property of tolerance proportionality [7] is asymptotically 
ensured through working for the local error with 
2 

(42) Est = - (4 - iT n A h>n ) r h (t n+1/2 ), A h>n — d Uh F h (t n , V h)n ) , 

where 7 is the stability coefficient of ROS3P. The common filter (4— , yT n Ah,n) serves to damp 
spurious stiff components which would otherwise be amplified through the i^-evaluations 
within r h (t n+1/2 ). 

Let D n = 1 1 -Est 1 1 and T ol n = T ol a + ToIr\\ Vh,n\\ with ToIa and ToIr given local tolerances. 
If D n > Tol n the step is rejected and redone. Otherwise the step is accepted and we advance 
in time. In both cases the new step size is determined by 

(43) r new = min(l.5,max(2/3,0.9r)) r n , r = (TolJ D^ 3 . 

After each step size change we adjust r new to r„ +1 = (T — t n )/|_(l + (T — t n )/r new )\ so as to 
guarantee to reach the end point T with a step of averaged normal length. The initial step 
size To is prescribed and is adjusted similarly. 

The linear error transport equations fl24|) and fl39]) are simultaneously solved by means of 
the implicit midpoint rule, which gives approximations th,n and ry^ n to the global time and 
spatial error at time t = t n . We use the implementations 

(h - 2T n Ah, n ) 5e n+ i = 2e h ^ n + §r n r(t n+ i/ 2 ) , 
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and 
(45) 



(4 - \T n A Kn ) 5r] n+1 = 2i] Kn - r n a h (t n+1/2 ) 



fj h>n+1 = 5l] n+1 - fj h>n . 

Clearly, the matrices Ah t n already computed within ROS3P can be reused. The spatial 
truncation error &2h(t) at t = t n +i/2 is given by 

(46) a 2h (t n+1 / 2 ) = ~ {R^h^h (Wl/2> Vh{tn+l/2)) ~ F 2h (Wl/2> tt^+l/j))) ■ 

Since Vh(t n +i/2) an d Fh(t n+ i/ 2 , Vh(t n+ i/ 2 ) are available from the computation of rh(t n+ i/ 2 ) in 
( |27|) . this requires only one function evaluation on the coarse grid. The vector a 2 h(t n+ i/ 2 ) on 
the coarse mesh is prolongated to the fine mesh and is then divided by 2 q = 4 if the neighboring 
fine grid points are equidistant, otherwise it is divided by 2 9_1 = 2. The remaining dh{t n+ i/ 2 ) 
on the fine mesh are computed by interpolation respecting the order of the neighboring 
spatial truncation errors. 

Due to freezing the coefficients in each time step, the second-order midpoint rule is a 
first-order method when interpreted for solving the linearised equations (I14p (or likewise 
the first variational equation) and fTTTj) . Thus if all is going well, we asymptotically have 
e h ,n = e h {t n ) + 0{r^ ax ) and fj h , n = Vh{t n ) + 0(T max h q max ) + 0(h q +^.). 

After computing the spatial truncation errors we can solve the discretized error transport 
equations fj45|) for all fjh, n - We shall distinguish between two different mesh adaptation 
approaches: (i) globally uniform and (ii) locally adaptive refinement. Although the uniform 
strategy may be less efficient, it is very easy to implement and therefore of special practical 
interest if software packages which do no allow dynamic adaptive mesh refinement are used. 

Uniform spatial refinement. Let Tol be a given tolerance. Then our aim is to guarantee 
||^/i(^)|| < Tol. From (1431) . we get an approximate value fj h>M for the spatial discretization 
error at T. If the desired accuracy is still not satisfied, i.e., ||%,m|| > Tol, we choose a new 
(uniform) spatial resolution 



(47) h 



new 




to account for achieving ||^h ne „(T)|| pa Tol. From h new we determine a new number of mesh 
points. The whole computation is redone with the new spatial mesh. 

Adaptive spatial refinement. The main idea of our local spatial mesh control is based 
on the observation that the principle of tolerance proportionality can be also applied to the 
spatial discretization error. Multiplying all <5/ l (t n+ i/ 2 ) in (j4"5l) by a certain constant multiplies 
all fjh, n +i by the same constant since fjh,o = 0. Set T ol^ = T oV\ + Tol^WV^nW where Tol% and 
Tol ft are given local tolerances and define a local estimator A n through 

(48) A\= £ 2h i \& l (t n+1/2 )\ 2 , 

i:Xi£F h 

where Fh denotes the set of all (fine) mesh points that do not belong to the coarse mesh. 
Remember we have second order of the spatial truncation error in these points. If A n < Tol" 
the mesh is only coarsened. Otherwise, if A n > Tol" the mesh is improved by refinement 
and coarsening as well. We set a to i = 0.9 Tol"/ \J~N and mark all Xi G Fh 

for refinement if v^«i(Wi/2) > OL toi 
and for coarsening if \fh~i oii{t ri+ i/ 2 ) < 0.1 a to i ■ 
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Grid adaptation is first performed for the coarse mesh and afterwards the fine mesh is 
constructed by halving each interval. If Xi is marked for refinement the corresponding coarse 
grid interval is halved. Grid points are only removed if there are two equidistant neighboring 
intervals the midpoints of which are marked for coarsening. Finally, the grid is smoothed 
such that 0.5 < hi/h^i < 2 everywhere. Data transfer from old to new meshes is done by 
cubic Hermite interpolation where the necessary first derivatives are determined from fourth 
order finite differences. 

After mesh adaptation the local time step is redone with the new mesh. The procedure 
is continued until first D n < Tol n and second A n < Tol% hold. The whole strategy aims at 
equidistributing the local values y/hi 0:1(^+1/2) • Asymptotically we get 

(50) ^J 2 £ al) =Uj2 °- 81( ^ <)2 ) - 0.9 Toll , 

\ i:Xi£F h / \ i:Xi&F h J 

where the factor 0.9 improves the robustness of the equidistribution principle. 

6. The control rules 

Like for the ODE case studied in [3] our aim is to provide global error estimates and to 
control the accuracy of the solution numerically computed to the imposed tolerance level. Let 
GToIa and GTol R be the global tolerances. Then we start with the local tolerances ToIa = 
GToIa, ToIr = GToIr, and in the spatially adaptive case also with Tol a A = C a GTolA, and 
Tol R = C a GTolji, where the factor C a > 1 ensures that the residual time error is small with 
respect to the spatial truncation error and therefore the use of fl38|) is justified. 

Suppose the numerical schemes have delivered an approximate solution Vh,M and global 
estimates eh,M and f\h,M for the time and spatial error at time tu — T. We then verify 
whether 

(51) \\^h,m\\ - CrCcontroiTolM, Tol M = GTol A + GTol R \\V hM \\, 

where C cont roi ~ 1, typically > 1, and Ct € (0, 1) denotes the fraction desired for the global 
time error with respect to the tolerance ToIm- If (EU) does not hold, the whole computation 
is redone over [0, T] with the same initial step Tq and the adjusted local tolerances 

(52) Tol A = Tol A - fac, Tol R = Tol R ■ fac, fac = C T Tol M /\\e h ,M\\- 

Based on tolerance proportionality, reducing the local error estimates with the factor fac will 
reduce e^(T) by fac [7]. 

If fl5T|) holds, we check whether 

(53) \\eh,M + Vh,M\\ < C con trolTolM- 

If it is true, the overall error E h {T) = Vh(T)—R h u(T) = eh(T)-\i] h (T) is considered small enough 
relative to the chosen tolerance and V^m is accepted. Otherwise, the whole computation is 
redone with the (already) adjusted tolerances ([52"]) and an improved spatial resolution. 

In the uniform case, we use the new mesh size computed from fj4T|) with Tol = (1 — 
Ct)ToIm- To check the convergence behaviour in space and therefore also the quality of 
the approximation of the spatial truncation error, we additionally compute the numerically 
observed order 

(54) g _ = log ( r ll^^/iog^). 

\\\Vh new ,M\\/ \hnewj 
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Step 


Control Algorithm with Uniform Refinement in Space 


Step 


Choose global tolerances GToIa and GTol R . 
Choose C T , C controh h , q, and r . 
Set local tolerances ToIa = GToIa and ToIr = GToIr. 
Set h = ho- 


Step 1 


Run numerical schemes to compute Vh,M-,^h,M^h,M- 
Compute ToIm = GToIa + GToIrWV^mW- 


Step 2 


IF ||e A) Af|| < C T C con troiTolM GOTO Step 3. 

ELSE set fac = CtToIm/W^mW, ToIa = ToIa • fac, ToIr = ToIr ■ fac 
and GOTO Step 1. 


Step 3 


IF \\e h)M + VkmW ^ CcontroiTol M GOTO Step 4. 

ELSE set h=j/(l- C T )Tol M /\\rjh,M\\ h and GOTO Step 1. 


Step 4 


IF h ^ h compute q num . 

ELSE set h = 2h, run numerical schemes again and compute then q num - 
IF q n um ~ q accept fine grid solution and STOP. 
ELSE set ho = 2h , h = h and GOTO Step 1. 



Table 1 . Algorithmic structure of the overall control strategy where uniform 
refinement in space is used. 



If qnum computed for the final run is not close to the expected value q used for our Richardson 
extrapolation, we reason that the approximation of the spatial truncation errors has failed 
due to a dominating global time error, which happens, e.g., if the initial spatial mesh is 
already too fine. Consequently, we coarsen the initial mesh by a factor two and start again. 
If the control approach stops without a mesh refinement, we perform an additional control 
run on the coarse mesh and compute q nU m from (|54j) with h new = 2/t. It turns out that 
this simple strategy works quite robustly, provided that the meshes used are able to resolve 
the basic behaviour of the solution. The algorithmic structure of our control strategy with 
uniform refinement in space is given in Table [TJ 

In the adaptive case, we choose new local tolerances 

(55) Toll = Tol a A ■ fac, Tol% = ToIr ■ fac, fac = (1 - C T )Tol M /\\fj h>M \\ , 

and the whole computation is redone over the interval [0, T\. Based on tolerance propor- 
tionality, reducing the local truncation error with the factor fac will reduce T)h(T) by fac. In 
Table [21 the algorithmic structure of our control strategy with adaptive refinement in space 
is displayed. Note that now the index h refers to a sequence of spatial meshes adapted at 
each time point t n . 

Summarizing, the first check flBT]) and the possibly second control computation serve to 
significantly reduce the global time error. This enables us to make use of the approximation 
( |38|) for the spatial truncation error, which otherwise could not be trusted. The second step 
based on suitable spatial mesh improvement attempts to bring the overall error down to the 
imposed tolerance. Using the sum of the approximate global time and spatial error inside 
the norm in (1531 . we take advantage of favourable effects of error cancellation. These two 
steps are successively repeated until the second check is successful. Additionally, if uniform 
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Step 


Control Algorithm with Adaptive Refinement in Space 


Step 


Choose global tolerances GToIa and GTol R . 
Choose C T , C 'control, C a , q, and r . 
Set local tolerances 

Tol a = GToIa, ToIr = GToIr, Tol A = C a GToIa, and Tol ft = C a GTol r. 
Choose initial spatial mesh. 


Step 1 


Run numerical schemes to compute V^m, e^M, ??h,M- 
Compute ToIm = GToIa + GToIr\\V^m\ • 


Step 2 


IF He^Afll < C T C contr oiTol M GOTO Step 3. 

ELSE set fac = CtToIm /\\&h m\\, ToIa = ToIa • fete, Tolp = Tola ■ fac 
and GOTO Step 1. 


Step 3 


IF \\eh,M + Vh,M\\ — CcontroiTolM accept solution and STOP. 
ELSE set 

fac=(l - C T )Tol M / 1| fjh,M ||, Tol a A = Tol a A ■ fac, Tol% = Tol R ■ fac 
and GOTO Step 1. 



Table 2. Algorithmic structure of the overall control strategy where adaptive 
refinement in space is used. 



mesh refinement is used we take into account the numerically observed order in space to 
assess the approximation of the spatial truncation error. 

7. Numerical illustrations 

To illustrate the performance of the global error estimators and the control strategy, we 
consider three test problems: (i) the highly stable heat equation with nonhomogeneous 
Neumann boundary conditions [1], (ii) the nonlinear convection-dominated Burgers' equa- 
tion [U [6], and (iii) the Allen-Calm equation modelling a diffusion-reaction problem [1]. 
Analytic solutions are known for all three problems. Uniform spatial refinement is studied 
for all three test cases. For the Allen-Cahn problem, these results are compared to those 
obtained with adaptive refinement. 

We set GTol A = GTol R = GTol for GTol = 10~ l , I = 2, . . . , 7 and start with one and 
the same initial step size r = 10~ 5 . Equally spaced meshes of 25, 51, 103, 207, 415, 831, 
and 1663 points are used as initial mesh. The control parameters introduced above for the 
control rules are CV = 1/3, C contro i = 1-2, and C a = 10. All runs were performed, but for 
convenience we only select a representative set of them for our presentation. 

We define the estimated global error E^m = &h,M + fjh,M at time t = T and set indicators 
O es t = || Eh,M ||/ 1 1-^/1(^)11 f° r the ratio of the estimated global error and the true global error, 
and Q ctr = ToIm l\\Eh(T)\\ for the ratio of the desired tolerance and the true global error. 
Thus, Q ctr > 1/ C C ontroi = 5/6 indicates control of the true global error. 

The tables of results contain the following quantities, Tol = ToIa = ToIr from (152|) . 
Tol a = Tol^ = Tol a R from fl35}, Tol M = GTol (1 + \\V KM \\) from flSTj), the estimated global 
error E^m, the estimated time error e^.M, and the estimated spatial truncation error fjh t M- 
Note that we always start with Tol = GTol in the first run. The ratios Q es t and Q c tr serve 
to illustrate the quality of the global error estimation and the control. If uniform refinement 
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in space is applied, the numerically observed order q num for the spatial error is given. It 
will be clear from the tables of results whether a tolerance-adapted run to control the global 
time error, a spatial mesh adaption step or an additional control run on a coarser grid was 
necessary. Especially, the latter is marked by a dashed line. 



Tol 


N 


Tol M 








©est 


@ctr 


Qnum 


1.00e-2 


25 


1.10e-2 


7.14e-4 


1.16e-4 


8.20e-4 


0.99 


15.27 




1.00e-2 


13 


1.10e-2 


3.27e-3 


1.24e-4 


3.38e-3 


0.99 


3.32 


2.04 


1.00e-3 


51 


1.10e-3 


1.68e-4 


1.97e-5 


1.86e-4 


1.00 


6.51 




1.00e-3 


25 


1.10e-3 


8.04e-4 


2.03e-5 


8.22e-4 


1.00 


1.36 


2.02 


1.00e-4 


103 


1.10e-4 


4.27e-5 


2.01e-6 


4.44e-5 


1.00 


2.57 




1.00e-4 


51 


1.10e-4 


1.85e-4 


1.96e-6 


1.86e-4 


1.00 


0.59 


2.01 


1.00e-5 


207 


1.10e-5 


1.07e-5 


1.89e-7 


1.08e-5 


1.00 


1.03 




1.00e-5 


103 


1.10e-5 


4.43e-5 H 


1.83e-7 H 


4.44e-5 


1.00 


0.25 


2.01 


1.00e-6 


415 


1.10e-6 


2.67e-6 


1.81e-8 


2.68e-6 


1.00 


0.41 




1.00e-6 


795 


1.10e-6 


7.14e-7 


1.83e-8 


7.28e-7 


1.00 


1.54 


2.00 


1.00e-7 


25 


1.10e-7 


8.24e-4 


1.24e-9 


8.24e-4 


1.00 


0.00 




1.00e-7 


2759 


1.10e-7 


5.91e-8 


1.60e-9 


6.03e-8 


1.00 


1.86 


2.01 


1.00e-7 


1663 


1.10e-7 


1.65e-7 


1.57e-9 


1.66e-7 


1.00 


0.67 




1.00e-7 


2505 


1.10e-7 


7.20e-8 


1.57e-9 


7.31e-8 


1.00 


1.53 


2.00 



Table 3. Selected data for the heat equation with Neumann boundary con- 
ditions. Uniform refinement in space is used. 



7.1. Heat equation with Neumann boundary conditions. This heat equation provides 
an example with inhomogeneous Neumann boundary conditions: 

(56) d t u = d xx u, 0<ar<1.0, 0<t<T = 0.2, 

and boundary conditions d x u = ir e -71 " 2 ' cos(7rx) at x = and x — 1. The initial condition is 
consistent with the analytic solution u(x,t) = e _7r 'sin(7ra). Although the solution is very 
stable, it is not easy to provide good error estimates as stated in [Q Ej. 

To approximate the inhomogeneous Neumann boundary conditions we introduce artificial 
mesh points x_i = —h and xn +2 = 1+h, discretize d x u(0) and d x u(l) by second order central 
differences, and use the approximate differential equation at the boundary to eliminate the 
artificial solution values. In consequence, we have second order in all mesh points and q = 2 
works also fine for interpolating the estimated spatial truncation error at the boundary. 

Due to the high stability of the problem the global time errors are much smaller than 
imposed local tolerances. So, control of the global time error is redundant here and control 
runs were only carried out in case of insufficient spatial resolutions. Table [3] shows results for 
various tolerances and initial meshes. The global error estimation and control appear to work 
very well for this problem, where the influence of the initial mesh points is less strong. This 
holds also for other combinations of tolerances and initial meshes. Note the high quality 
of the estimator E^m (and therefore also of %,m), showing that the derivative boundary 
condition is well resolved within the Richardson extrapolation. For the runs with tolerances 
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GTol = 10 2 , 10 3 , 10 4 , 10 5 , the order of the spatial convergence was successfully checked 
with a second run on the coarse mesh, that is, we can trust the first run. 



Tol 


N 


Tol M 




IpftjAf || 


\\Vh, M || 


©est 


®ctr 


Qnum 


1.00e-2 


51 


1.93e-2 


4.30e-3 


1.86e-3 


2.86e-3 


1.08 


4.87 




1.00e-2 


25 


1.93e-2 


1.29e-2 


2.21e-3 


1.14e-2 


0.99 


1.48 


2.00 


1.00e-3 


51 


1.93e-3 


2.83e-3 


1.54e-4 


2.74e-3 


0.99 


0.68 




1.00e-3 


75 


1.93e-3 


1.36e-3 


1.48e-4 


1.28e-3 


1.00 


1.42 


2.00 


1.00e-4 


51 


1.93e-4 


2.73e-3 


1.09e-5 


2.73e-3 


0.98 


0.07 




1.00e-4 


239 


1.94e-4 


1.32e-4 


1.05e-5 


1.27e-4 


1.00 


1.46 


2.00 


1.00e-5 


51 


1.93e-5 


2.73e-3 


1.08e-6 


2.73e-3 


0.98 


0.01 




1.00e-5 


757 


1.94e-5 


1.32e-5 


1.02e-6 


1.27e-5 


1.00 


1.47 


2.00 


1.00e-6 


51 


1.93e-6 


2.73e-3 


1.08e-7 


2.73e-3 


0.98 


0.00 




1.00e-6 


2391 


1.94e-6 


1.32e-6 


9.29e-8 


1.28e-6 


1.00 


1.47 


2.00 


1.00e-7 


51 


1.93e-7 


2.73e-3 


1.10e-8 


2.73e-3 


0.98 


0.00 




1.00e-7 


7563 


1.94e-7 


1.31e-7 


8.57e-9 


1.28e-7 


1.00 


1.47 


2.00 



Table 4. Selected data for Burgers' equation with 51 initial mesh points. 
Uniform refinement in space is used. 



7.2. Burgers' equation. The second problem is the nonlinear Burgers' equation 

(57) d t u = e d xx u - ud x u , 0<x<1.0, 0<t<T=1.0, 

where e = 0.015 is used in the experiments. Dirichlet boundary conditions and initial 
conditions are consistent with the analytic solution defined by 

f - a s f A n + 5r 2 + 10r 3 

(58) U(X,t) = ; r 

y J K ! 10(n + r 2 + r 3 ) 

where r x {x) = e 0A5x / £ , r 2 {t,x) = e - 01 ( 10 + 6 *+ 25 *)/ £ , and r 3 (t) = e - 025 ( 6 - 5 + 9 - 9i )/ £ . 

In Table S] we present results for all tolerances used and the 51-point initial mesh. The use 
of a relatively coarse mesh at the beginning is the natural choice in practice. No adaptation in 
time is necessary, which is mainly due to the small first time step and the maximum factor 
1.5 which is allowed in fT43|) for a step size enlargement. For the tolerance GTol = 10~ 2 , 
the numerical solution is accepted since the corresponding control run shows g n Mm~2, the 
expected value. Remarkably excellent estimates are obtained for higher tolerances. Here, 
control is always achieved after one spatial mesh improvement. 

7.3. The Allen-Cahn equation. The third problem is the bi-stable Allen-Cahn equation 
which is defined by 

(59) d t u = 10~ 2 d xx u + lOOu (I - u 2 ) , < x < 2.5 , < t < T = 0.5 , 

with the initial function and Dirichlet boundary values taken from the exact wave front 
solution u(x,t) = (1 + gM^-"*))-^ \ — 50 y/2, a = 1.5 y/2. This problem was also used in 

HI- 

First we apply uniform refinement in space. Table [5] reveals a high quality of the global 
error estimation and also the control process works quite well. Let us pick one exemplary run 
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Tol 


N 


T 1 Ctl n /T 
± ULM 


II ^h,M || 


II °h,M || 


Will 71 t\\ 

II >lh,M || 


w est 


8 * 

^ctr 


n 

Hnum 


l.UUe-Z 


iuo 


Z.UOc-Z 


1 84e> n 


1 K„ 1 


1 Q8o n 
i .yoe-u 


Q 8Q 

y.oy 


n 1 1 

U. 1 1 




4 fiQo 4 


IUO 


9 D^P 9 


R 7Sp_1 

O. f Oc- 1 


1 9Rp 3 

1 .ZUc-O 


^ 7Qp 1 
o. / ye-i 


z.uy 


n i n 






u i i 


9 D9p 9 

Z.UZc _ Z 


fi 04p 3 


1 lip Q 
1 . 1 le-O 


7 1 3 


1 1 Q 
i . i y 


^ Q8 
o.yo 


9 "\A 


l.UUe-Z 


41 ^ 


9 f19c 9 

Z.UZc-Z 


7 RQp 9 
i . uye-z 


1 44p 1 


fi 7^0 9 

U. 1 <JC-Z 


o.uo 


n 8n 

U.oU 




zL fifip A 


41 ^ 


9 D9p 9 

Z.UZc _ Z 


1 Sfip 9 


1 lip Q 


1 Q7p 9 
i . y i e-z 


1 9^ 
1 .zo 


1 ^4 

1 .O^l 






907 
ZU ( 


9 fl^o 9 
Z.UOe-Z 


y . 1 1 e-z 


1 . lOe-O 


Q 9Qo 9 
y .zye-z 


1 47 
1 i 


n ^9 

u.oz 


9 94 


l.UUe-O 


9D7 
ZU l 


Z.Uoe-O 


Q 89o 9 

y.oze-z 


9 Q7o 
z.y i e-o 


1 .U ie-1 


1 fin 

l.DU 


u.uo 




9 97o /I 

z.z < e-^t 


9D7 
ZU ( 


Z.Uoe-O 


o.oue-z 


A Q^o /I 


8 8^o 9 
O.oOe-Z 


1 

i.oy 


u.uo 




9 97c. 4 

z.z i 


1 fiS'? 
IDoO 


9 D9o 3 
Z.UZe-O 


fi 1 /lo A 


4 71 o 4 


1 flQp 1 

i .uye-o 


1 1 1 


^ fi7 
0.0 ( 


9 i n 

Z. 1U 


l.UUe-o 


SQ1 

ool 


z.Uze-o 


z.zoe-o 


z.o / e-o 


K 1 On Q 

o.l ze-o 


1 QQ 
l.OO 


1 1 n 
l.ly 




2.35e-4 


831 


2.02e-3 


4.01e-3 


4.91e-4 


4.50e-3 


1.12 


0.57 




2.35e-4 


1521 


2.02e-3 


8.42e-4 


4.90e-4 


1.33e-3 


1.12 


2.68 


2.02 


1.00e-4 


1663 


2.02e-4 


8.89e-4 


1.86e-4 


1.08e-3 


1.07 


0.24 




3.63e-5 


1663 


2.02e-4 


9.88e-4 


6.14e-5 


1.05e-3 


1.05 


0.21 




3.63e-5 


4643 


2.02e-4 


7.30e-5 


6.14e-5 


1.34e-4 


1.04 


2.89 


2.00 



Table 5. Selected data for the Allen-Calm problem. Uniform refinement in 
space is used. 



Tol 


Tol a 


N M 


Tol M 


Eh,m\ 


Zh,Af\ 


Vh,M 


©est 


©ctr 


1.00e-2 
4.81e-4 


1.00e-l 
1.00e-l 


245 
247 


2.01e-2 
2.01e-2 


1.05e-l 
8.54e-3 


1.39e-l 
1.13e-3 


3.42e-2 
9.67e-3 


3.21 
1.26 


0.61 

2.98 


1.00e-3 
2.35e-4 


1.00e-2 
1.00e-2 


483 
481 


2.01e-3 
2.01e-3 


1.26e-3 
9.72e-4 


2.86e-3 
4.84e-4 


1.59e-3 
1.46e-3 


1.26 
1.11 


2.01 
2.30 


1.00e-4 
3.62e-5 


1.00e-3 
1.00e-3 


1839 
1839 


2.01e-4 
2.01e-4 


9.08e-5 
5.49e-5 


1.85e-4 
6.06e-5 


9.45e-5 
1.16e-4 


1.63 
0.92 


3.62 
3.36 


1.00e-4 
3.62e-5 
3.62e-5 


1.00e-2 
1.00e-2 
9.68e-4 


481 
483 
1839 


2.01e-4 
2.01e-4 
2.01e-4 


1.23e-3 
1.32e-3 
5.48e-5 


1.85e-4 
6.09e-5 
6.07e-5 


1.41e-3 
1.38e-3 
1.15e-4 


1.08 
1.06 
0.92 


0.18 
0.16 
3.36 


1.00e-4 
3.62e-5 
3.62e-5 


1.00e-l 
1.00e-l 
1.55e-3 


243 
243 
1809 


2.01e-4 
2.01e-4 
2.01e-4 


8.66e-3 
8.55e-3 
5.68e-5 


1.84e-4 
6.08e-5 
6.06e-5 


8.84e-3 
8.61e-3 
1.17e-4 


1.15 
1.12 
0.92 


0.03 
0.03 
3.25 



Table 6. Selected data for the Allen-Calm problem. Adaptive refinement in 
space is used. 



out to explain the overall control strategy in more detail. Starting with GTol = Tol = 10~ 3 
and 831 mesh points, which corresponds to the fourth simulation, the numerical scheme 
delivers global error estimates H^mII = 2.87 x 10~ 3 and ||%,m|| = 5.12 x 10~ 3 for the time 
and spatial error of the approximate solution V^m at the final time tu = T. The first check 
for the time error estimate H^mH < CTC CO ntro{TolM = 8.08 x 10" 4 fails and we adjust the 
local tolerances by a factor fac = CtToIm/W^hmW = 2-35 x 10 -1 , which yields the new 
Tol = 2.35 10~ 4 . The computation is then redone. Due to the tolerance proportionality, in 
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the second run the time error is significantly reduced and the inequality H^mII < 8.08 x 10 4 
is now valid. We proceed with checking H-Et^mH < C contro iTolu = 2.42 x 10~ 3 , which is still 
not true. From (jSJ), we compute a new number of spatial mesh points N = 1521. Finally, 
the third run is successful and with the numerically observed spatial order q num = 2.02 the 
numerical solution is accepted. 

The ratios for Q est = HE^^/H/H^^)!! lie between 1.04 and 1.23, after the control runs. 
Control of the global error, that is H-E^T)!! < C contro iTolM, is in general achieved after two 
steps (one step to adjust the time grid and one step to control the space discretization), 
whereas the efficiency index Q ctr = ToIm l\\Eh{T)\\ is close to three. This results from a 
systematic cancellation effect between the global time and spatial error, which is not taken 
into account when computing h new from (T4T1) . 

Next we consider locally adaptive spatial grid enhancement instead of globally uniform 
adaptation. Within each time step the grid is adapted by refinement and coarsening, based on 
an eqidistribution principle, until A n < Tol% = Tol a (l + \\Vh in \\) holds. This yields a sequence 
of non-uniform meshes. Let Nm denote the number of adaptive grids points obtained at the 
final time T. The first three runs in Table O correspond to our standard setting C a = 10, 
i.e., Tol a = 10 Tol. In this case, after adjusting the local tolerances for the time integration 
no further run with higher tolerances in space is necessary. To demonstrate the robustness 
of the algorithm, we select two additional runs with C a = 10 l ,l = 2,3, for GTol = 10 -4 . 
In both cases, coarser meshes are used at the beginning and a second control run has to be 
done to decrease the spatial discretization error. The resulting adaptive spatial meshes are 
comparable. Control of the global error is always achieved. The estimation process works 
again quite well. Compared to the uniform case, significantly less spatial degrees of freedoms 
are needed to reach the desired tolerances. The reduction rate varies between 40% and 70%. 



8. Summary 

We have developed an error control strategy for finite difference solutions of parabolic 
equations, involving both temporal and spatial discretization errors. The global time error 
strategy discussed in [4] appears to provide an excellent starting point for the development 
of such an algorithm. The classical ODE approach based on the first variational equation 
and the principle of tolerance proportionality is combined with an efficient estimation of 
the spatial error and mesh adaptation to control the overall global error. Two approaches 
have been presented to handle spatial mesh improvement: (i) globally uniform refinement 
and (ii) local refinement and coarsening based on an equidistribution principle. Inspired 
by [I], we have used Richardson extrapolation to approximate the spatial truncation error 
within the method of lines. Our control strategy aims at balancing the spatial and temporal 
discretization error in order to achieve an accuracy imposed by the user. 

The key ingredients are: (i) linearised error transport equations equipped with sufficiently 
accurate defects to approximate the global time error and global spatial error and (ii) uniform 
or adaptive mesh refinement and local error control in time based on tolerance proportionality 
to achieve global error control. For illustration of the performance and effectiveness of our 
approach, we have implemented second-order finite differences in one space dimension and 
the example integrator ROS3P [3j. On the basis of three different test problems we could 
observe that our approach is very reliable, both with respect to estimation and control. 



Global Error Estimation and Control 



15 



9. Acknowledgement 

This work has been supported in part by the German Research Foundation (DFG) by the 
grant SFB568. 

References 

[1] M. Berzins (1988), Global error estimation in the method of lines for parabolic equations, SIAM J. Sci. 

Stat. Comput. 9, pp. 687-703. 
[2] J. Lang (2000), Adaptive Multilevel Solution of Nonlinear Parabolic PDE Systems. Theory, Algorithm 

and Applications, Lecture Notes in Computational Science and Engineering, Vol. 16, Springer. 
[3] J. Lang, J.G. Vcrwcr (2001), ROS3P - An accurate third-order Rosenbrock solver designed for parabolic 

problems, BIT 41, pp. 731-738. 
[4] J. Lang, J.G. Verwer (2007), On global error estimation and control for initial value problems, SIAM 

J. Sci. Comput. 29, pp. 1460-1475. 
[5] S. Larsson, V. Thomee (2005), Partial Differential Equations with Numerical Methods, Texts in Applied 

Mathematics, Vol. 45, 2nd printing, Springer. 
[6] L. Lawson, M. Berzins, P.M. Dew (1991), Balancing space and time errors in the method of lines for 

parabolic equations, SIAM J. Sci. Stat. Comput. 12, pp. 573-594. 
[7] L.F. Shampinc (1994), Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New 

York. 

[8] R.D. Skeel (1986), Thirteen ways to estimate global error, Numer. Math. 48, pp. 1-20. 
[9] J.W. Thomas (1995), Numerical Partial Differential Equations. Finite Difference Methods, Text in 
Applied Mathematics 22, Springer. 

Kristian Debrabant, Department of Mathematics, Technische Universitat Darmstadt, 
Dolivostr. 15, D-64293 Darmstadt, Germany 

E-mail address: debrabant@mathematik.tu-darmstadt.de 

Jens Lang, Department of Mathematics, Technische Universitat Darmstadt, Dolivostr. 15, 
D-64293 Darmstadt, Germany 

E-mail address: lang@mathematik.tu-darmstadt.de 



